
import numpy as np
import obspy as ob

from obspy import Trace
from scipy import signal
import tqdm
from obspy.core import UTCDateTime

from obspy.signal.invsim import cosine_taper, waterlevel

def get_tapered_traces(data:np.ndarray, dt:float, L_Taper=1, WIN_TYPE='cosine')-> np.ndarray:
    """
    对数据的最后一个维度上进行灭尖操作
    data:np.array([-1,t])
    dt     : 采样间隔
    L_Taper : 灭尖长度 , 单位为s
    WIN_TYPE: 窗函数类型 ('hann', 'cosine', 'tukey'等)

    return : tapered_traces
    """
    shape_raw = list(data.shape)
    nt = shape_raw[-1]
    data_2d = data.reshape([-1,nt]).copy()
    ne,_ = data_2d.shape
    
    taper_length = int(L_Taper/dt)
    taper_window = signal.get_window(WIN_TYPE, taper_length*2)
    # 前向taper
    data_2d[:,:taper_length] *= taper_window[:taper_length]
    # 反向taper (避免两端突变)
    data_2d[:,-taper_length:] *= taper_window[taper_length:]

    return data_2d.reshape(shape_raw)

def get_tapered_slices(data:np.ndarray, t:np.ndarray, T_win:list, L_Taper=1, WIN_TYPE='hann')-> np.ndarray:
    """
    对多维数据进行截取操作，叠加灭尖操作
    
    该函数从输入的多维数据中提取指定时间窗口内的数据，并应用窗函数进行边缘平滑处理，
    以减少截断效应。主要用于地震数据处理中的时间窗口截取。
    
    Parameters:
    data : numpy.ndarray
        输入的多维数据，最后一个维度为时间轴，形状为(..., nt)
    t : numpy.ndarray
        时间轴数组，形状为(nt,)，与data的最后一个维度对应
    T_win : list
        时间窗口范围 [t_start, t_end]，单位为秒
    L_Taper : float, optional
        灭尖长度，单位为秒，默认为1秒
        如果小于等于0，则不进行灭尖操作
    WIN_TYPE : str, optional
        窗函数类型，默认为'hann'
        支持的类型包括'hann', 'cosine', 'tukey'等
    
    Returns:
    sliced_data : numpy.ndarray
        截取并经过灭尖处理的数据，形状为(..., ntc)
        其中ntc为截取后的时间点数
    tc : numpy.ndarray
        截取后的时间轴数组，形状为(ntc,)
    
    Examples:
    >>> import numpy as np
    >>> data = np.random.random((3, 1000))  # 3个台站，1000个时间点
    >>> t = np.linspace(0, 10, 1000)  # 0-10秒的时间轴
    >>> sliced_data, tc = get_tapered_slices(data, t, [2, 8], L_Taper=0.5)
    >>> print(sliced_data.shape)  # (3, 600) 假设采样率为60Hz
    >>> print(tc[0], tc[-1])  # 2.0, 8.0
    """

    shape_raw = list(data.shape)
    nt = shape_raw[-1]
    data_2d = data.reshape([-1,nt])
    ne,_ = data_2d.shape
    data_cut = []

    dt = t[1]-t[0]
    sp = int((T_win[0]-t[0])/dt)
    ep = int((T_win[1]-t[0])/dt)
    if sp<0 or ep>nt:
        raise ValueError(f'时间窗口超出数据范围,{T_win=}')

    data_cut = data_2d[:,sp:ep].copy()

    tc = t[sp:ep]
    if L_Taper>0:
        data_cut = get_tapered_traces(data_cut, dt, L_Taper, WIN_TYPE=WIN_TYPE)

    data_cut = np.array(data_cut).reshape(shape_raw[:-1]+[len(tc)])

    return data_cut, tc

def get_tapered_slices_multiWIN(data:np.ndarray, t:np.ndarray, T_win:list, L_Taper=1, WIN_TYPE='hann')-> np.ndarray:
    """
    对多维数据进行截取操作，对不同的T时间窗口进行不同窗函数的设置
    
    该函数从输入的多维数据中提取指定时间窗口内的数据，并应用窗函数进行边缘平滑处理，
    以减少截断效应。主要用于地震数据处理中的时间窗口截取。
    
    Parameters:
    data : numpy.ndarray
        输入的多维数据，最后一个维度为时间轴，形状为(..., nt)
    t : numpy.ndarray
        时间轴数组，形状为(nt,)，与data的最后一个维度对应
    T_win : list
        时间窗口范围 [[t_start, t_end]]，单位为秒
    L_Taper : float, optional
        灭尖长度，单位为秒，默认为1秒
        如果小于等于0，则不进行灭尖操作
    WIN_TYPE : str, optional
        窗函数类型，默认为'hann'
        支持的类型包括'hann', 'cosine', 'tukey'等
    
    Returns:
    sliced_data : numpy.ndarray
        截取并经过灭尖处理的数据，形状为(..., ntc)
        其中ntc为截取后的时间点数
    tc : numpy.ndarray
        截取后的时间轴数组，形状为(ntc,)
    
    Examples:
    >>> import numpy as np
    >>> data = np.random.random((3, 1000))  # 3个台站，1000个时间点
    >>> t = np.linspace(0, 10, 1000)  # 0-10秒的时间轴
    >>> sliced_data, tc = get_tapered_slices(data, t, [2, 8], L_Taper=0.5)
    >>> print(sliced_data.shape)  # (3, 600) 假设采样率为60Hz
    >>> print(tc[0], tc[-1])  # 2.0, 8.0
    """

    shape_raw = list(data.shape)
    nt = shape_raw[-1]
    data_2d = data.reshape([-1,nt])
    ne,_ = data_2d.shape
    data_cut = []

    dt = t[1]-t[0]
    sp = int((T_win[0]-t[0])/dt)
    ep = int((T_win[1]-t[0])/dt)
    if sp<0 or ep>nt:
        raise ValueError(f'时间窗口超出数据范围,{T_win=}')

    data_cut = data_2d[:,sp:ep].copy()

    tc = t[sp:ep]
    if L_Taper>0:
        data_cut = get_tapered_traces(data_cut, dt, L_Taper, WIN_TYPE=WIN_TYPE)

    data_cut = np.array(data_cut).reshape(shape_raw[:-1]+[len(tc)])

    return data_cut, tc

def get_tapered_slice(trace:np.ndarray, t:np.ndarray, T_win:list, L_Taper:list)-> np.ndarray:
    """
    使用 ObsPy 的 taper 和 trim 方法对数据进行截取操作
    return : sliced_trace, times
    """
    # 将数据转为 ObsPy Trace 对象
    tr = Trace(data=np.asarray(trace))
    tr.stats.sampling_rate = 1.0 / (t[1] - t[0])  # 计算采样率
    tr.stats.starttime = t[0]
    # Step 1: 截取时间窗口（自动处理边界）
    trimmed_tr = tr.trim(starttime= tr.stats.starttime+T_win[0]-t[0], 
                        endtime=tr.stats.starttime+T_win[1]-t[0], 
                        pad=True, fill_value=0)
    # Step 2: 对截取后的数据应用灭尖
    trimmed_tr.taper(max_percentage=None, max_length=L_Taper[0], side="left")
    trimmed_tr.taper(max_percentage=None, max_length=L_Taper[1], side="right")

    return trimmed_tr.data, trimmed_tr.times()+T_win[0]

def array2Trace(data, dt, t0=UTCDateTime()) -> ob.Trace:
    '''
    将多维数据转为 ObsPy Trace 对象
    return : ob.Trace
    '''
    tr = Trace(data=np.asarray(data))
    tr.stats.sampling_rate = 1.0 / dt  # 计算采样率
    tr.stats.starttime = t0

    return tr
def get_traces_spec(traces, dt,F_LIM=[2,4]):

    '''
    对多维数据最后一个维度进行快速频谱计算
    traces: np.array([-1,nt])
    dt: 采样间隔
    F_LIM: 频谱范围 [fs,fe]
    return : f_freq, amp
    '''
    nt = traces.shape[-1]
    traces = traces.reshape([-1,nt])

    fs,fe = F_LIM
    spec1 = np.fft.fft(traces, axis=-1)
    f_freq = np.fft.fftfreq(nt, dt)
    
    phi = np.angle(spec1)
    amp = np.abs(spec1)
    
    fsp = np.argmin(np.abs(f_freq-fs))
    fep = np.argmin(np.abs(f_freq-fe))
    
    f_freq = f_freq[fsp:fep]
    amp = amp[:,fsp:fep]

    return f_freq,amp

def safe_decimate(data, original_sr, target_sr, taper_percentage=0.01, taper_type='hann'):
    """
    安全降采样函数（带边界taper和抗混叠滤波）
    
    参数：
    - data: data
    - original_sr: 原始采样率 (Hz)
    - target_sr: 目标采样率 (Hz)
    - taper_percentage: 两端taper的百分比 (0.05 = 5%)
    - taper_type: taper类型 ('hann', 'cosine', 'tukey'等)
    
    返回：
    - decimated_trace: 降采样后的Trace对象
    """
    trace = Trace(data=np.asarray(data))
    trace.stats.sampling_rate = original_sr  # 计算采样率
    # 检查采样率合法性
    original_sr = trace.stats.sampling_rate
    if target_sr >= original_sr:
        raise ValueError("目标采样率必须低于原始采样率")
    
    # 计算降采样因子
    decimation_factor = int(original_sr / target_sr)
    if not np.isclose(original_sr / decimation_factor, target_sr, rtol=1e-6):
        raise ValueError("采样率需为整数倍关系")

    # 应用边界taper（避免振铃效应）
    trace.taper(max_percentage=taper_percentage, type=taper_type)

    # 分阶段降采样（当降采样因子过大时）
    current_sr = original_sr
    while current_sr > target_sr:
        # 计算当前阶段的最大安全降采样因子
        stage_factor = min(decimation_factor, 10)  # 单次最多10倍
        
        # 计算实际目标采样率
        stage_sr = current_sr / stage_factor
        
        # 执行降采样（包含内置抗混叠滤波）
        trace.decimate(factor=stage_factor, no_filter=False)
        
        # 更新当前采样率
        current_sr = stage_sr
        decimation_factor = int(current_sr / target_sr)

    return trace

def fast_decimate(data, t_raw, target_dt, L_Taper=1, taper_type='hann', fmin=0.3):
    '''
    对多维数据最后一个维度进行快速降采样（带边界taper）
    
    参数：
    - data: data, [-1,nt]
    - t_raw: time axis
    - target_dt: 目标采样间隔（s）
    - L_Taper: taper的长度（s）
    - taper_type: taper类型 ('hann', 'cosine', 'tukey'等)
    
    返回：
    - decimated_data, t: 降采样后的Trace对象
    '''
    shape_raw = list(data.shape)
    nt = shape_raw[-1]
    data_2d = data.reshape([-1,nt])

    dt_raw = t_raw[1]-t_raw[0]
    original_sr = 1/dt_raw
    target_sr = 1/target_dt
    
    # taper
    data_2d = get_tapered_traces(data_2d.copy(), dt_raw, L_Taper=L_Taper, WIN_TYPE=taper_type)

    # 低通滤波
    fs = target_sr/2
    data_2d = filtfilt(data_2d, dt_raw, [fmin, fs*0.95], N_CORE=None, type='bandpass')
    
    # 降采样
    decimation_factor = int(original_sr / target_sr)
    if not np.isclose(original_sr / decimation_factor, target_sr, rtol=1e-6):
        raise ValueError("采样率需为整数倍关系")
    data_2d = data_2d[:,::decimation_factor]
    t_n = t_raw[::decimation_factor]

    # 恢复原始尺寸
    data_nd = data_2d.reshape(shape_raw[:-1]+[-1])
    

    return data_nd,t_n

def filtfilt(data:np.ndarray, dt, FS, N_CORE=None, type='bandpass', BAR=False,order=1):
    '''
    对多维地震数据进行滤波
    data: 地震数据 (numpy数组)
    dt: 采样间隔 (s)
    FS: 滤波频率 (Hz)
    N_CORE: 使用的核数 (int)
    type: 滤波类型 ('bandpass', 'lowpass', 'highpass')
    
    return data: filtered_data
    '''
    f_max = 1/dt/2
    if type=='bandpass':
        fs,fe = FS
        b,a = signal.butter(order, [fs/f_max,fe/f_max],"bandpass")
    elif type=='lowpass':
        fe = FS
        b,a = signal.butter(1, fe/f_max,"low")
    elif type=='highpass':
        fs = FS
        b,a = signal.butter(order, fs/f_max,"high")
    else:
        raise ValueError(f'No such type value {type}')
    shape_raw = data.shape
    nt = shape_raw[-1]
    data = data.copy()
    data_2d = data.reshape([-1,nt])
    
    ne,_ = data_2d.shape

    if N_CORE is None:
        # iter = range(ne)
        # if BAR:
        #     iter=tqdm.tqdm(iter)
        # for ie in iter:
        #     trace =  data_2d[ie,:]
        data_2d = signal.filtfilt(b,a,data_2d, axis=-1) 
                                  
        return data_2d.reshape(shape_raw)
    else:
        from multiprocessing.pool import Pool

        N_job=N_CORE*2
        NE_per_core = ne//N_job+1
        input_args = []
        for j in range(N_job):
            s_j = j*NE_per_core
            e_j = min(ne,(j+1)*NE_per_core)
            data_j = data_2d[s_j:e_j,:]
            input_args.append([b,a,data_j])    

        results=[]
        with Pool(N_CORE) as pool:
            results =results + pool.starmap(__filtfilt_parallel_task,input_args)

        data_2d=np.zeros([ne,nt], dtype='float32')
        for j in range(N_job):
            data_j = results[j]
            s_j = j*NE_per_core
            e_j = min(ne,(j+1)*NE_per_core)
            data_2d[s_j:e_j,:] = data_j
        print('multiprocesses end')
        del results

        return data_2d.reshape(shape_raw)
def __filtfilt_parallel_task(b,a,data):
    ne,nt = data.shape
    data_new = np.zeros_like(data)
    for i in range(ne):
        data_new[i,:] = signal.filtfilt(b,a,data[i,:])
    return data_new

def safe_filter(mytrace, dt, ftype='bandpass', zerophase=True,  **paras):
    mytrace = Trace(data=np.asarray(mytrace))
    mytrace.stats.sampling_rate = 1/dt

    mytrace.filter(ftype, **paras, zerophase=zerophase)
    return mytrace.data
def get_mean_trace(data_1d, dt, L):
    '''
    对多维地震数据进行均值滤波
    data_1d: 地震数据 (numpy数组)
    dt: 采样间隔 (s)
    L: 滤波长度 (s)
    
    return data: filtered_data
    '''

    FL = int(L/dt)//2*2
    mean_filter = np.ones([FL]) 
    npts = len(data_1d)
    d_mean = np.zeros_like(data_1d)
    d_mean[(FL//2):(npts-FL//2+1)] = np.convolve(data_1d, mean_filter,'valid')/len(mean_filter)

    return d_mean

# ... existing code ...
def get_traces_attributes(data, t, WIN_L, WIN_stride=1, ATTR = 'mean'):
    '''
    对多维地震数据,获取分段最大值、平均振幅等能量信息
    data: 地震数据 (numpy数组), 最后一个维度是时间轴
    t: data所对应的时间
    L: 长度 (s)
    stride: 步长，(s)
    
    return data: filtered_data (仅保持有效值)
    '''
    
    # 保持原始形状
    original_shape = list(data.shape)
    nt = len(t)
    dt = t[1]-t[0]
    # 将数据重塑为2D数组，最后一个维度是时间轴
    data_2d = data.reshape([-1, nt])
    ne, _ = data_2d.shape
    
    WLP = int(WIN_L/dt)
    SWP = int(WIN_stride/dt)
    
    # 使用for循环逐行处理，避免strided_app可能的内存问题
    nrows = (nt - WLP) // SWP+1
    # 初始化结果数组
    result = np.zeros((ne, nrows, WLP))
    tt = t[np.arange(0, nt-WLP, SWP)]+SWP/2*dt
    # 逐行处理每条地震道数据
    for i in range(ne):
        for j in range(nrows):
            start_idx = j * SWP
            end_idx = start_idx + WLP
            result[i, j, :] = data_2d[i, start_idx:end_idx]
    
    if ATTR=='mean':
        d_mean_2d = np.mean(np.abs(result), axis=2)
        d_mean = d_mean_2d.reshape(original_shape[:-1]+[-1])
        return d_mean, tt
    elif ATTR=='max':
        d_max_2d = np.max(np.abs(result), axis=2)
        d_max = d_max_2d.reshape(original_shape[:-1]+[-1])
        return d_max, tt
    elif ATTR=='E':
        d_E_2d = np.sqrt(np.mean(result**2, axis=2))
        d_E = d_E_2d.reshape(original_shape[:-1]+[-1])
        return d_E, tt
    else:
        raise ValueError(f'Invalid attribute {ATTR}, not in [mean, max, E]')

def get_traces_sqrtmean(data, t, WIN_L):
    '''
    对多维地震数据,获取分段最大值、平均振幅等能量信息
    data: 地震数据 (numpy数组), 最后一个维度是时间轴
    t: data所对应的时间
    WIN_L: 长度 (s)
    
    return data: filtered_data (仅保持有效值)
    '''
    
    # 保持原始形状
    original_shape = list(data.shape)
    nt = len(t)
    dt = t[1]-t[0]
    # 将数据重塑为2D数组，最后一个维度是时间轴
    data_2d = data.reshape([-1, nt])
    ne, _ = data_2d.shape
    
    WLP = int(WIN_L/dt)
    N = nt//WLP
    tN = np.reshape(t[:N*WLP],[N, WLP])
    tN = np.mean(tN, axis=-1)

    data_2d_N = np.reshape(data_2d[:,:N*WLP],[-1,N, WLP])
    
    d_E_2d = np.sqrt(np.mean(data_2d_N**2, axis=2))*np.sqrt(2)
    
    d_E = d_E_2d.reshape(original_shape[:-1]+[-1])
    return d_E, tN
    

def stack_traces_in_bins(traces, x, bin_size=32, x_min=None, x_max=None):
    """
    将迹线按指定距离范围和箱大小进行分箱堆叠
    
    Parameters:
    traces: array, 输入的迹线数据
    x: array, 对应的距离位置
    bin_size: int, 箱的大小(默认32米)
    x_min: int, 距离范围最小值(默认-1000米)
    x_max: int, 距离范围最大值(默认1000米)
    
    Returns:
    tuple: (stacked_traces, stacked_x, stacked_names)
    """
    print(f'Stacking traces in {bin_size}m bins...')
    x_min =  x.min() if x_min is None else x_min
    x_max =  x.max() if x_max is None else x_max
    bin_edges = np.arange(x_min, x_max + bin_size, bin_size)
    bin_centers = bin_edges[:-1] + bin_size/2

    # 初始化新的trace数组
    stacked_traces = []
    stacked_x = []
    stacked_names = []

    for i in range(len(bin_edges) - 1):
        left_edge = bin_edges[i]
        right_edge = bin_edges[i+1]
        
        # 找到在当前bin范围内的trace
        indices = np.where((x >= left_edge) & (x < right_edge))[0]
        
        if len(indices) > 0:
            # 将这些trace叠加（求平均）
            stacked_trace = np.mean(traces[indices, :], axis=0)
            stacked_traces.append(stacked_trace)
            stacked_x.append(bin_centers[i])
            stacked_names.append(f"Bin_{left_edge:.0f}_{right_edge:.0f}_n{len(indices)}")

    # 转换为numpy数组
    stacked_traces = np.array(stacked_traces)
    stacked_x = np.array(stacked_x)
    
    return stacked_traces, stacked_x, stacked_names


def spectral_whitening(data, dt, fs,fe):
    '''
    进行谱白化
    data: 地震数据 (numpy数组)
    dt: 采样间隔 (s)
    fs: 滤波频率范围 (Hz)
    fe: 滤波频率范围 (Hz)
    '''
    pass
    
def whiten_msnoise(data, dt, freqmin, freqmax, FL=1):
    """This function takes 1-dimensional *data* timeseries array,
        goes to frequency domain using fft, whitens the amplitude of the spectrum
        in frequency domain between *freqmin* and *freqmax*
        and returns the whitened fft.
        :type data: :class:`numpy.ndarray`
        :param data: Contains the 1D time series to whiten
        :type delta: float
        :param delta: The sampling frequency of the `data`
        :type freqmin: float
        :param freqmin: The lower frequency bound
        :type freqmax: float
        :param freqmax: The upper frequency bound
        :FL: float
        :param FL: The length of corner window, HZ

        :rtype: :class:`numpy.ndarray`
        :returns: The FFT of the input trace, whitened between the frequency bounds

        from https://github.com/fablindner/glseis/blob/master/msnoise_move2obspy.py
    """

    Nfft = len(data)
    freqVec = np.fft.rfftfreq(Nfft, d=dt)
    df = freqVec[1]-freqVec[0]
    Napod = int(FL/df)


    J = np.where((freqVec >= freqmin) & (freqVec <= freqmax))[0]
    low = J[0] - Napod
    if low <= 0:
        low = 1

    porte1 = J[0]
    porte2 = J[-1]
    high = J[-1] + Napod
    #if high > Nfft / 2:
    #    high = int(Nfft // 2)

    FFTRawSign = np.fft.rfft(data, Nfft)

    # Left tapering:
    FFTRawSign[0:low] *= 0
    FFTRawSign[low:porte1] = np.cos(
        np.linspace(np.pi / 2., np.pi, porte1 - low)) ** 2 * np.exp(
        1j * np.angle(FFTRawSign[low:porte1]))
    # Pass band:
    FFTRawSign[porte1:porte2] = np.exp(1j * np.angle(FFTRawSign[porte1:porte2]))
    # Right tapering:
    FFTRawSign[porte2:high] = np.cos(
        np.linspace(0., np.pi / 2., high - porte2)) ** 2 * np.exp(
        1j * np.angle(FFTRawSign[porte2:high]))
    #FFTRawSign[high:Nfft + 1] *= 0
    FFTRawSign[high:] *= 0

    wdata = np.real(np.fft.irfft(FFTRawSign, Nfft))

    return wdata

def moving_ave(A,N):
    '''
    this Numba compiled function does running smooth average for an array.
    PARAMETERS:
    ---------------------
    A: 1-D array of data to be smoothed
    N: integer, it defines the half window length to smooth

    RETURNS:
    ---------------------
    B: 1-D array with smoothed data
    '''
    A = np.concatenate((A[:N],A,A[-N:]),axis=0)
    B = np.zeros(A.shape,A.dtype)

    tmp=0.
    for pos in range(N,A.size-N):
        # do summing only once
        if pos==N:
            for i in range(-N,N+1):
                tmp+=A[pos+i]
        else:
            tmp=tmp-A[pos-N-1]+A[pos+N]
        B[pos]=tmp/(2*N+1)
        if B[pos]==0:
            B[pos]=1
    return B[N:-N]

def whiten_noisepy(data, dt,fs,fe,smooth_N, freq_norm,Napod=10):
    '''
    This function takes 1-dimensional timeseries array, transforms to frequency domain using fft,
    whitens the amplitude of the spectrum in frequency domain between *freqmin* and *freqmax*
    and returns the whitened fft.
    PARAMETERS:
    ----------------------
    data: numpy.ndarray contains the 1D time series to whiten
    fft_para: dict containing all fft_cc parameters such as
        dt: The sampling space of the `data`
        freqmin: The lower frequency bound
        freqmax: The upper frequency bound
        smooth_N: integer, it defines the half window length to smooth
        freq_norm: whitening method between 'one-bit' and 'RMA'
    RETURNS:
    ----------------------
    FFTRawSign: numpy.ndarray contains the FFT of the whitened input trace between the frequency bounds
    '''

    import scipy
    # load parameters
    delta   = dt
    freqmin = fs
    freqmax = fe
    smooth_N  = smooth_N
    freq_norm = freq_norm

    # Speed up FFT by padding to optimal size for FFTPACK
    if data.ndim == 1:
        axis = 0
    elif data.ndim == 2:
        axis = 1

    Nfft = data.shape[-1]

    Nfft = int(Nfft)
    freqVec = scipy.fftpack.fftfreq(Nfft, d=delta)[:Nfft // 2]
    J = np.where((freqVec >= freqmin) & (freqVec <= freqmax))[0]
    low = J[0] - Napod
    if low <= 0:
        low = 1

    left = J[0]
    right = J[-1]
    high = J[-1] + Napod
    if high > Nfft/2:
        high = int(Nfft//2)

    FFTRawSign = scipy.fftpack.fft(data, Nfft,axis=axis)
    # Left tapering:
    if axis == 1:
        FFTRawSign[:,0:low] *= 0
        FFTRawSign[:,low:left] = np.cos(
            np.linspace(np.pi / 2., np.pi, left - low)) ** 2 * np.exp(
            1j * np.angle(FFTRawSign[:,low:left]))
        # Pass band:
        if freq_norm == 'phase_only':
            FFTRawSign[:,left:right] = np.exp(1j * np.angle(FFTRawSign[:,left:right]))
        elif freq_norm == 'rma':
            for ii in range(data.shape[0]):
                tave = moving_ave(np.abs(FFTRawSign[ii,left:right]),smooth_N)
                FFTRawSign[ii,left:right] = FFTRawSign[ii,left:right]/tave
        # Right tapering:
        FFTRawSign[:,right:high] = np.cos(
            np.linspace(0., np.pi / 2., high - right)) ** 2 * np.exp(
            1j * np.angle(FFTRawSign[:,right:high]))
        FFTRawSign[:,high:Nfft//2] *= 0

        # Hermitian symmetry (because the input is real)
        FFTRawSign[:,-(Nfft//2)+1:] = np.flip(np.conj(FFTRawSign[:,1:(Nfft//2)]),axis=axis)
    else:
        FFTRawSign[0:low] *= 0
        FFTRawSign[low:left] = np.cos(
            np.linspace(np.pi / 2., np.pi, left - low)) ** 2 * np.exp(
            1j * np.angle(FFTRawSign[low:left]))
        # Pass band:
        if freq_norm == 'phase_only':
            FFTRawSign[left:right] = np.exp(1j * np.angle(FFTRawSign[left:right]))
        elif freq_norm == 'rma':
            tave = moving_ave(np.abs(FFTRawSign[left:right]),smooth_N)
            FFTRawSign[left:right] = FFTRawSign[left:right]/tave
        # Right tapering:
        FFTRawSign[right:high] = np.cos(
            np.linspace(0., np.pi / 2., high - right)) ** 2 * np.exp(
            1j * np.angle(FFTRawSign[right:high]))
        FFTRawSign[high:Nfft//2] *= 0

        # Hermitian symmetry (because the input is real)
        FFTRawSign[-(Nfft//2)+1:] = FFTRawSign[1:(Nfft//2)].conjugate()[::-1]

    RawSign = np.real(scipy.fftpack.ifft(FFTRawSign,Nfft,axis=axis))
    
    return RawSign